function PV = func_polynomial_vector(nv,np)

% nv = number of variables {1,2,3,4}
% np = degree of polynomials

C = zeros((np+1)^nv,nv);

%% Numbe of variables = 1
if nv == 1
    for j1=0:np
        jj = j1;
        C(jj,:) = j1;
    end
    PV = C;
end
        
%% Number of variables = 2
if nv == 2
    for j2=0:np
        for j1=0:np
            jj = j1+1+(np+1)*j2;
            C(jj,:) = [j1,j2];
        end
    end
    sumC = sum(C,2);
    I = sumC < np+1;
    PV = C(I,:);
end

%% Number of variables = 3
if nv == 3
    for j3=0:np
        for j2=0:np
            for j1=0:np
                jj = j1+1 + (np+1)*j2 + (np+1)^2*j3;
                C(jj,:) = [j1,j2,j3];
            end
        end
    end
    sumC = sum(C,2);
    I = sumC < np+1;
    PV = C(I,:);
end

%% Number of variables = 4
if nv == 4
    for j4=0:np
        for j3=0:np
            for j2=0:np
                for j1=0:np
                    jj = j1+1 + (np+1)*j2 + (np+1)^2*j3 + (np+1)^3*j4;
                    C(jj,:) = [j1,j2,j3,j4];  % all the possible combinations of 0,1,2,3
                end
            end
        end
    end
    sumC = sum(C,2);     % column vector containing sums of each row
    I    = sumC < np+1;  % indices for sum of each column < np+1
    PV   = C(I,:);       % all the possible combinations of 0,1,2,3, that are less than np+1
end

